Matching
Data
We’ll load up the Lalonde data. Our goal here is to see if we can replicate or at least get close to the correct estimate in Lalonde’s 1986 paper using observational data. According to his model, the actual effect for the job training program was a $1,794 increase in the subject’s wages in 1978.
As usual a good starting point here is to get some basic descriptive stats on our variables.
lalonde|>
select(treat, educ, married, nodegree, race)|>
mutate(treat = factor(treat, labels=c("control", "treatment")),
married = factor(married, labels=c('single','married')),
nodegree = factor(nodegree, labels=c('degree','no degree'))
)|>
ggpairs(aes(fill=treat, color=treat), upper='blank') +
theme_bw() +
scale_fill_brewer(palette='Dark2') +
scale_color_brewer(palette='Dark2')lalonde|>
select(treat, age, re74, re75, re78)|>
mutate(treat = factor(treat, labels=c("control", "treatment")))|>
ggpairs(aes(fill=treat, color=treat), upper='blank') +
theme_bw() +
scale_fill_brewer(palette='Dark2') +
scale_color_brewer(palette='Dark2')Regression
We’ll start with regression based estimates. The first model includes no controls, and the second includes controls for age, education, marital status, past income, degree status, and race:
# the uncontrolled regression
mod0<-lm(re78 ~ treat, data=lalonde)
mod1<-lm(re78 ~ treat + age + educ + married + re74 + re75+ nodegree +race,data=lalonde)
modelsummary(list("treatment only" =mod0,
"treatment + controls" = mod1),
estimate = "{estimate}",
fmt = fmt_significant(),
statistic ='conf.int',
conf_level = .95,
gof_omit = 'F|RMSE|R2$|AIC|Log.Lik.', # remove some model stats
# add a title
title = "DV: Wages"
)| treatment only | treatment + controls | |
|---|---|---|
| (Intercept) | 6984 | -1174 |
| [6276, 7693] | [-5998, 3650] | |
| treat | -635 | 1548.2 |
| [-1926, 655] | [13.9, 3082.6] | |
| age | 13.0 | |
| [-50.8, 76.8] | ||
| educ | 403.9 | |
| [91.9, 716.0] | ||
| married | 407 | |
| [-959, 1772] | ||
| re74 | 0.2964 | |
| [0.1819, 0.4108] | ||
| re75 | 0.2315 | |
| [0.0261, 0.4370] | ||
| nodegree | 260 | |
| [-1404, 1924] | ||
| racehispan | 1740 | |
| [-261, 3740] | ||
| racewhite | 1241 | |
| [-269, 2750] | ||
| Num.Obs. | 614 | 614 |
| R2 Adj. | 0.000 | 0.135 |
| BIC | 12712.0 | 12666.1 |
Matching methods
Mahalanobis Scores
You can get Mahalanobis score matches by using method="nearest" and distance='mahalanobis'
ma_matched <- matchit(treat ~ age + educ + married + race +
nodegree + re74 + re75,
data = lalonde,
method = "nearest",
distance= 'mahalanobis',
replace=T,
ratio=5
)Once we have our matchit object, the Cobalt package gives us some options for balance checking. After examining the results, you might want to go back and re-specify the model.
Balance Measures
Type Diff.Un Diff.Adj M.Threshold
age Contin. -0.3094 0.2904 Not Balanced, >0.1
educ Contin. 0.0550 -0.0403 Balanced, <0.1
married Binary -0.3236 0.0314 Balanced, <0.1
race_black Binary 0.6404 0.0086 Balanced, <0.1
race_hispan Binary -0.0827 0.0000 Balanced, <0.1
race_white Binary -0.5577 -0.0086 Balanced, <0.1
nodegree Binary 0.1114 0.0130 Balanced, <0.1
re74 Contin. -0.7211 0.0817 Balanced, <0.1
re75 Contin. -0.2903 0.1488 Not Balanced, >0.1
Balance tally for mean differences
count
Balanced, <0.1 7
Not Balanced, >0.1 2
Variable with the greatest mean difference
Variable Diff.Adj M.Threshold
age 0.2904 Not Balanced, >0.1
Sample sizes
Control Treated
All 429. 185
Matched (ESS) 74.59 185
Matched (Unweighted) 196. 185
Unmatched 233. 0
bal.plot(ma_matched)love.plot(ma_matched)Access the matched data with the match_data function:
ma_data<-match_data(ma_matched)
head(ma_data)Note that a fairly large number of observations have been discarded here:
nrow(ma_data)[1] 381
We can use the results in our regression model. Be sure to include the weights!
ma_fit<-lm(re78 ~ treat + age + educ + married + race +nodegree + re74 + re75,
data = ma_data,
weight = weights )Finally, we want to estimate our effects with robust standard errors:
avg_comparisons(ma_fit,
variables = "treat",
vcov = 'HC3'
)
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
1152 900 1.28 0.201 2.3 -612 2915
Term: treat
Type: response
Comparison: 1 - 0
modelsummary(list(
"Mahalanobis" = ma_fit
),
estimate = "{estimate}",
fmt = fmt_significant(),
statistic ='conf.int',
conf_level = .95,
vcov ='HC3',
gof_omit = 'F|RMSE|R2$|AIC|Log.Lik.', # remove some model stats
# add a title
title = "DV: Wages"
)| Mahalanobis | |
|---|---|
| (Intercept) | -2890 |
| [-10778, 4998] | |
| treat | 1152 |
| [-618, 2921] | |
| age | -20.8 |
| [-117.8, 76.2] | |
| educ | 703 |
| [117, 1290] | |
| married | 160 |
| [-1937, 2257] | |
| racehispan | 1718 |
| [-1243, 4678] | |
| racewhite | 1246 |
| [-588, 3080] | |
| nodegree | 885 |
| [-2036, 3807] | |
| re74 | 0.0772 |
| [-0.2979, 0.4524] | |
| re75 | 0.200 |
| [-0.192, 0.593] | |
| Num.Obs. | 381 |
| R2 Adj. | 0.035 |
| BIC | 7982.9 |
| Std.Errors | HC3 |
Coarsened Exact matching
For CEM, we’ll often want to set manually set some of “cutpoints”. There’s no set rules here. Visual inspection of the predictors can help, and the MatchIt package has some defaults that can be sensible, but this is partly a judgement call: what values of the independent variables seem like they should be grouped together?
By default, observations are matched exactly on factor variables. However, you can override this by setting passing a list of levels to the grouping argument:
cem_match <- matchit(treat ~ age + educ + married + race +
nodegree + re74 + re75, data = lalonde,
method = "cem",
estimand ='ATE',
cutpoints =list(educ = c(0,9, 12, 14)),
grouping = list(race = list(c("white", "hispan"),
c("black")))
)
summary(cem_match)
Call:
matchit(formula = treat ~ age + educ + married + race + nodegree +
re74 + re75, data = lalonde, method = "cem", estimand = "ATE",
cutpoints = list(educ = c(0, 9, 12, 14)), grouping = list(race = list(c("white",
"hispan"), c("black"))))
Summary of Balance for All Data:
Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
age 25.8162 28.0303 -0.2419 0.4400 0.0813
educ 10.3459 10.2354 0.0448 0.4959 0.0347
married 0.1892 0.5128 -0.7208 . 0.3236
raceblack 0.8432 0.2028 1.6708 . 0.6404
racehispan 0.0595 0.1422 -0.2774 . 0.0827
racewhite 0.0973 0.6550 -1.4080 . 0.5577
nodegree 0.7081 0.5967 0.2355 . 0.1114
re74 2095.5737 5619.2365 -0.5958 0.5181 0.2248
re75 1532.0553 2466.4844 -0.2870 0.9563 0.1342
eCDF Max
age 0.1577
educ 0.1114
married 0.3236
raceblack 0.6404
racehispan 0.0827
racewhite 0.5577
nodegree 0.1114
re74 0.4470
re75 0.2876
Summary of Balance for Matched Data:
Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
age 21.1665 20.7539 0.0451 0.8970 0.0144
educ 10.1082 10.1923 -0.0340 0.7899 0.0153
married 0.0619 0.0619 -0.0000 . 0.0000
raceblack 0.5876 0.5876 0.0000 . 0.0000
racehispan 0.1314 0.0940 0.1256 . 0.0375
racewhite 0.2809 0.3184 -0.0946 . 0.0375
nodegree 0.7423 0.7423 0.0000 . 0.0000
re74 472.0782 742.1934 -0.0457 1.2640 0.0478
re75 466.2790 677.2163 -0.0648 1.1178 0.0549
eCDF Max Std. Pair Dist.
age 0.1294 0.1133
educ 0.1186 0.3455
married 0.0000 0.0000
raceblack 0.0000 0.0000
racehispan 0.0375 0.2567
racewhite 0.0375 0.1933
nodegree 0.0000 0.0000
re74 0.2828 0.0847
re75 0.2699 0.1654
Sample Sizes:
Control Treated
All 429. 185.
Matched (ESS) 86.99 32.26
Matched 112. 82.
Unmatched 317. 103.
Discarded 0. 0.
Now we can plot our results to check out how well the matching performed:
plot(cem_match,
type = "density",
interactive = TRUE,
which.xs = ~age +
married + re75)Balance Measures
Type Diff.Un Diff.Adj M.Threshold
age Contin. -0.2419 0.0451 Balanced, <0.1
educ Contin. 0.0448 -0.0340 Balanced, <0.1
married Binary -0.3236 -0.0000 Balanced, <0.1
race_black Binary 0.6404 0.0000 Balanced, <0.1
race_hispan Binary -0.0827 0.0375 Balanced, <0.1
race_white Binary -0.5577 -0.0375 Balanced, <0.1
nodegree Binary 0.1114 0.0000 Balanced, <0.1
re74 Contin. -0.5958 -0.0457 Balanced, <0.1
re75 Contin. -0.2870 -0.0648 Balanced, <0.1
Balance tally for mean differences
count
Balanced, <0.1 9
Not Balanced, >0.1 0
Variable with the greatest mean difference
Variable Diff.Adj M.Threshold
re75 -0.0648 Balanced, <0.1
Sample sizes
Control Treated
All 429. 185.
Matched (ESS) 86.99 32.26
Matched (Unweighted) 112. 82.
Unmatched 317. 103.
bal.plot(cem_match)love.plot(cem_match)Finally, we’ll want to use the match_data function to extract the matched data and then estimate our model with the weights included:
cem_data<-match_data(cem_match)
cem_fit <- lm(re78 ~
treat + age + educ + married + re74 +re75+
nodegree+ race , data = cem_data, weights = weights)avg_comparisons(cem_fit, variables = "treat", vcov ='HC3', cluster=~subclass)
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
1672 1158 1.44 0.149 2.7 -597 3941
Term: treat
Type: response
Comparison: 1 - 0
modelsummary(list(
"Mahalanobis" = ma_fit,
"Coarsened" = cem_fit
),
estimate = "{estimate}",
fmt = fmt_significant(),
statistic ='conf.int',
conf_level = .95,
vcov ='HC3',
gof_omit = 'F|RMSE|R2$|AIC|Log.Lik.', # remove some model stats
# add a title
title = "DV: Wages"
)| Mahalanobis | Coarsened | |
|---|---|---|
| (Intercept) | -2890 | -4693 |
| [-10778, 4998] | [-16885, 7499] | |
| treat | 1152 | 1672 |
| [-618, 2921] | [-612, 3956] | |
| age | -20.8 | 106 |
| [-117.8, 76.2] | [-104, 317] | |
| educ | 703 | 608 |
| [117, 1290] | [-164, 1379] | |
| married | 160 | -2533 |
| [-1937, 2257] | [-6329, 1263] | |
| racehispan | 1718 | 957 |
| [-1243, 4678] | [-1046, 2959] | |
| racewhite | 1246 | 2391 |
| [-588, 3080] | [-853, 5636] | |
| nodegree | 885 | 826 |
| [-2036, 3807] | [-2735, 4387] | |
| re74 | 0.0772 | -0.0964 |
| [-0.2979, 0.4524] | [-1.1358, 0.9430] | |
| re75 | 0.200 | 0.687 |
| [-0.192, 0.593] | [-0.214, 1.588] | |
| Num.Obs. | 381 | 194 |
| R2 Adj. | 0.035 | 0.061 |
| BIC | 7982.9 | 3999.8 |
| Std.Errors | HC3 | HC3 |
Inverse probability weighting
Finally, we can try using inverse probability weighting on the propensity scores instead of matching. Here’s how we might do that using mostly base R functions:
# Step 1: Estimate propensity scores
fit1 <- glm(treat ~age + educ + nodegree +
married + race + re74 + re75 ,
family = binomial, data = lalonde)
lalonde$propensity <- predict(fit1, type = "response")
# calculate the inverse
lalonde<-lalonde|>
mutate(ipw = (treat/ propensity) + ((1 - treat) / (1 - propensity)))
# Step 2: Fit weighted outcome model
m <- lm(re78 ~ treat , data = lalonde, weight = ipw)
summary(m)
Call:
lm(formula = re78 ~ treat, data = lalonde, weights = ipw)
Weighted Residuals:
Min 1Q Median 3Q Max
-42083 -6606 -2284 4979 77674
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6422.8 397.4 16.161 <2e-16 ***
treat 224.7 577.7 0.389 0.697
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9864 on 612 degrees of freedom
Multiple R-squared: 0.0002471, Adjusted R-squared: -0.001386
F-statistic: 0.1513 on 1 and 612 DF, p-value: 0.6975
avg_comparisons(m, variables = "treat", wts = lalonde$ipw, vcov = 'HC0')
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
225 909 0.247 0.805 0.3 -1558 2007
Term: treat
Type: response
Comparison: 1 - 0
Or we can do this using the WeightIt package:
W <- weightit(treat ~ age + educ + nodegree +
married + race + re74 + re75,
data = lalonde, method = "glm",
estimand = "ATE")
fit <- lm_weightit(re78 ~ treat, data = lalonde,
weightit = W)
summary(fit, ci = TRUE)
Call:
lm_weightit(formula = re78 ~ treat, data = lalonde, weightit = W)
Coefficients:
Estimate Std. Error z value Pr(>|z|) 2.5 % 97.5 %
(Intercept) 6422.8 353.4 18.177 <1e-06 5730.3 7115.4 ***
treat 224.7 876.2 0.256 0.798 -1492.6 1942.0
Standard error: HC0 robust (adjusted for estimation of weights)
Or trying using a different weighting method:
W <- weightit(treat ~ age + educ + nodegree +
married + race + re74 + re75,
data = lalonde, method = "ebal",
estimand = "ATE")
fit <- lm_weightit(re78 ~ treat, data = lalonde,
weightit = W)
summary(fit, ci = TRUE)
Call:
lm_weightit(formula = re78 ~ treat, data = lalonde, weightit = W)
Coefficients:
Estimate Std. Error z value Pr(>|z|) 2.5 % 97.5 %
(Intercept) 6327.5 326.2 19.400 <1e-06 5688.3 6966.8 ***
treat 951.7 1232.6 0.772 0.44 -1464.2 3367.5
Standard error: HC0 robust (adjusted for estimation of weights)
avg_comparisons(fit, variables = "treat")
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
952 1233 0.772 0.44 1.2 -1464 3368
Term: treat
Type: probs
Comparison: 1 - 0